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We use non-equilibrium molecular dynamics simulations to study the adverse role of defects in- 
cluding isotopic impurities on the thermal conductivity of carbon nanotubes, graphene and graphene 
nanoribbons. We find that even in structurally perfect nanotubes and graphene, isotopic impurities 
reduce thermal conductivity by up to one half by decreasing the phonon mean free path. An even 
larger thermal conductivity reduction, with the same physical origin, occurs in presence of structural 
defects including vacancies and edges in narrow graphene nanoribbons. Our calculations reconcile 
results of former studies, which differed by up to an order of magnitude, by identifying limitations 
of various computational approaches. 
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I. INTRODUCTION 

With increasing performance of microprocessors, ris- 
ing heat evolution poses a serious problem. 1 To prevent 
damage, excess heat is conducted away to a heat sink 
using interconnects with high thermal conductivity. In 
diamond, which is used for this purpose and which con- 
ducts heat by phonons, isotopic impurities reduce its ex- 
cellent thermal conduction by up to one halfP^ The ini- 
tial prediction that thermal conductivity of perfect car- 
bon nanotubes and graphene monolayers (not graphite) 
should be similar and even surpass the diamond values 4 
was subsequently confirmed experimentally, albeit with a 
large scatter in the observed valuesPThe added benefit of 
nanotubes and graphene is their dual role as thermal con- 
ductors and active elements in electronic circuits. Unlike 
in heavier elements, the 13 C/ 12 C mass ratio does modify 
the phonon spectra of graphitic nanostructures signifi- 
cantly, causing a large reduction in thermal conductivity 
of systems with both isotopesP^ 

To obtain microscopic understanding of factors lim- 
iting thermal conductivity in graphitic nanostructures, 
we perform large-scale non-equilibrium molecular dy- 
namics (MD) simulations of defective carbon nanotubes, 
graphene and graphene nanoribbons. We determine 
the temperature-dependent thermal conductivity of 12 C- 
based systems as a function of 13 C concentration and 
compare the effect of isotopic impurities to that of diva- 
cancies. We show that, depending on temperature, the 
thermal conductivity of 13 C x 12 Ci- x nanostructures may 
be quenched by up to one half with respect to isotopically 
pure systems. Even at low concentrations, atomic vacan- 
cies quench thermal conductance more efficiently than 
isotopic impurities. Whereas freely suspended graphene 
monolayers conduct heat almost as well as isolated car- 
bon nanotubes, edge scattering reduces significantly the 
thermal conductivity of graphene nanoribbons. Our cal- 
culations reconcile results of former studies, which dif- 
fered by up to an order of magnitude, by identifying lim- 
itations of various computational approaches. 

Since in carbon nanostructures the electronic density 



of states at the Fermi level is either zero (diamond, 
graphene) or very small (nanotubes, graphene nanorib- 
bons), thermal transport in these systems is dominated 
by phonons. According to Fourier law, the thermal con- 
ductivity A is given by the heat current dQ/dt through 
area A in response to a temperature gradient dT/dz as 



^dQ 
A~dt 



dT 

dz 
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The phonon component of the thermal conductivity, 
which is dominant, is the product A = (l/3)cyv s (l), 
where cy is the specific heat per volume, v s the speed 
of sound, and (I) is the phonon mean free path. Rigid 
interatomic bonds in both sp 2 and sp 3 carbon struc- 
tures are responsible for a very high speed of sound 
v s and hard phonon modes, which translate into a 
high Debye frequency and large value of the specific 
heat cy. In isotopically pure monocrystalline diamond 
and carbon nanotubes, the phonon mean free path (I) 
may approach a large fraction of a micrometer, giving 
rise to record thermal conductivity values^ as large as 
^40, 000 Wm^K" 1 near T^lOO K. 

Since presence of defects, including isotopic impuri- 
ties and atomic vacancies of different types, can not be 
avoided in realistic systems, it is imperative to under- 
stand their role in thermal conductivity. We expect de- 
fects to play only a minor role in changing the speed 
of sound and specific heat, but to reduce drastically the 
phonon mean free path and thus the value of A. We be- 
lieve that the large scatter in the observed datsP comes 
not only from the extreme difficulty to measure this quan- 
tity in excellent thermal conductors, but more impor- 
tantly due to different types and concentrations of de- 
fects in different samples. Since controlling defects on 
the nanometer scale is nearly impossible experimentally, 
computer simulations provide a welcome alternative to 
understand the effect of particular defects on thermal 
conductivity. 
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II. METHOD 

To simulate computationally the conduction of heat, 
we make use of large-scale non-equilibrium molecular dy- 
namics (NEMD) simulations which had been used 
successfully to predict thermal conductivity of nanotubes 
and graphene. 4 The alternative way to calculate A us- 
ing direct MD simulations based on Eq. ([I]) requires ap- 
plying a thermostat that maintains a finite temperature 
difference AT across a finite distance Az. To prevent 
artifacts, Az must be larger than the phonon mean free 
path of up to 1 /im, which is computationally impractica- 
ble. A second alternative, which does not suffer from this 
limitation,^ is based on the Green-Kubo formula^ that 
relates A to the time-averaged autocorrelation function of 
the heat flux in the system. As shown earlier, 4 this time 
average converges very slowly in an equilibrium MD sim- 
ulation and depends sensitively on the initial conditions, 
making extensive ensemble averaging a necessary re- 
quirement that is computationally extremely demanding 
for systems of interest hereP Due to these problems,^ 
thermal conductivity calculations of graphitic nanostruc- 
tures based on direct MD,^ NEMDp> the Green- 
Kubo formalism] 14 * 18 * 1 ^ or the L andauer non-equilibrium 
Green's function formalis m 1 20 * 21 * have arrived at inconsis- 
tent results that differed by up to an order of magnitude 
and thus need to be revisited. 

Our computational approach 8 combines the Green- 
Kubo form ula^ with non-equilibrium molecular 
dynamic d 9 * 10 * in a computationally efficient manner P*l 
The dynamics of the system is driven by forces 



dt 2 



-VtU - (rriiVi + AFi 



(2) 



which act on individual atoms i. Here, r\ is the position 
and Vi the velocity of atom i with mass m^. 

The first term is the gradient of the total potential en- 
ergy U of the system, representing the force caused by 
interatomic interactions. To reduce the unusually high 
computational requirements, we represent U by the Ter- 
soff bond-order potential,^ which reproduces well the op- 
timum structure as well as vibration spectra of graphene 
and related nanostructures! 2 -^ This potential also for- 
mally allows a decomposition of the total potential energy 
into potential energies Ui of individual atoms, U = J2 i 
which will be of use in the following. 

The s econd term describes a Nose-Hoover 
thermostat! 25 * 26 * that represents coupling of the sys- 
tem to a heat bath at temperature T. The dynamics of 
the generalized coordinate ( of the heat bath is governed 

by 



d( 1 / 3jV-6 



(3) 



where Ek is the kinetic energy of the TV-atom system 
and Q the thermal inertia of the heat bath. 

The third term is a small fictitious force that acts as 
a perturbation, driving the system out of equilibrium to 



generate a heat flux, and is given by 

AF< = Ae,F e - f«(ry -F e ) + ^ E E { ^(v jk -Y e ) • 



(4) 

Here, F e is a vector parameter (with the dimension of 
inverse length) representing the strength of the perturba- 
tion. Tij = Tj — Ti and Aei = — (e) is the excess energy 
of atom z, where = m^ 2 /2 + U{ is its instantaneous 
energy and (e) = l/Af^e^. = —ViUj represents the 
contribution to the force on atom i stemming from its 
interaction with atom j, where is the gradient with 
respect to the position of atom i. 

The heat flux in the system is then given by 

d N 



(5) 

Setting F e = F e z, the ^-component of the heat flux 
can be obtained using the simplified expressions^ 

1 



J 2 



V,;AF,; 



(6) 



which leads to the thermal conductivity A along the z- 
direction 



A = lim lim 

F e ^0*^OO F e TV 



(7) 



This approach to determine A has been shown to 
be equivalent to that obtained using the Green- 
Kubo formulapmi yet is computationally much less 
demanding! 2 ^ 

To determine A in defect-free and defective nanotubes, 
graphene and graphene nanoribbons, we integrated the 
equations of motion using At = 0.2 fs as time step. We 
used the fifth-order predictor-corrector algorithm 2 ^ for 
Eq. ([2| and the fourth-order algorithm to integrate the 
coupled Eq. We used Q = 10 a.m.u.-A 2 for the ther- 
mal inertial of the thermostat, which allowed for efficient 
thermalization while not disturbing significantly the dy- 
namics of the system. The number of time steps needed 
for a reliable time average of the heat flux in Eq. ^ var- 
ied depending on the system, the temperature, and the 
value of F e . We found that only 500,000 time steps were 
sufficient to reach convergence for F e > 10- 3 A" 1 , but 
for smaller values 10" 5 A -1 < F e < 10" 3 A -1 we used 
up to 2xl0 6 time steps covering a 0.4 ns time period. 
The estimated 10-20% error in the extrapolation of our 
results towards F e — ?>0 is shown by the error bars of A in 
Figs. fTpl 



III. RESULTS AND DISCUSSION 

A. Thermal conductivity of carbon nanotubes 

Our results for the thermal conductivity of perfect 
and defective (10, 10) carbon nanotubes are presented 
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FIG. 1. (Color online) Thermal conductivity A of perfect and defective (10, 10) carbon nanotubes as a function of temperature 
T. (a) A in structurally perfect 13 C x 12 Ci- x nanotubes with varying isotopic composition in comparison to pure 12 C nanotubes 
with 0.5% missing atoms forming divacancy defects, (b) Details of (a) on a reduced temperature scale, (c) A of 12 C nanotubes 
with divacancies, presented in (a) and (b), on an expanded A scale, (d) Depiction of the nanotube unit cell containing different 
types of defects. The lines in (a)-(c) are guides to the eye. 



in Fig. [T] We used periodic boundary conditions with 
a large unit cell containing 400 C atoms, depicted in 
Fig. [jjd). Focussing first on structurally perfect nan- 
otubes with the isotopic composition 13 C x 12 Ci_ x , we 
present in Figs. [TJa) and [TJb) the thermal conductivity 
of isotopically pure ( 12 C and 13 C) nanotubes and isotope 
mixtures with x = 0.05 and x = 0.10 in the temperature 
range up to 600 K. In all systems, the low-temperature 
behavior of A is dominated by that of cy, which, along 
with A, approaches zero for T— ^0 K and then gradually 
increases with increasing temperature. After reaching its 
maximum, which occurs near T^lOO K in nanotubes, 
A decreases again due to the decreasing phonon mean 
free path, caused by increasing structural disorder at 
high temperatures. We find the thermal conductivity to 
be highest in isotopically pure 12 C nanotubes, with iso- 
topically pure 13 C nanotubes reaching almost the same 
value. We represented nanotubes with an isotopic mix- 
ture x = 0.05 and x = 0.10 by randomly distributing 
13 C atoms across the 12 C lattice. Our results for these 
mixtures indicate that for T < 300 K the thermal con- 
ductivity may decrease by up to 30% with respect to 
the value in isotopically pure lattices due to the strong 
reduction of the phonon mean free pathP^ Close inspec- 
tion of our results reveals that thermal conductivity of 
nanotubes remains almost unaffected by the presence of 
13 C isotopic impurities at temperatures T>300 K, where 
phonon-phonon scattering seems to dominate the mean 
free path reduction. 

To find out the relative importance of isotopic impuri- 
ties and structural defects, we also studied thermal con- 
ductivity in isotopically pure 12 C (10, 10) nanotubes con- 
taining a small fraction of atomic vacancies. We focussed 



on divac ancie s, which are more stable than monatomic 
vacancies, 1 29 1 30 1 and considered one single divacancy per 
400-atom unit cell. Our results for A in this system, 
shown in Figs. [TJa)-[TJc), indicate that even a very low 
concentration of structural defects may quench thermal 
conductivity by roughly an order of magnitude. These re- 
sults support our intuition that vacancies scatter phonons 
very efficiently and reduce the phonon mean free path 
even more than isotopic impurities. Whereas the reduc- 
tion of A by 97% at T^lOO K is extremely large, the 
relative role of structural defects decreases at higher tem- 
peratures, reaching a value of 83% at 300 K and 67% at 
600 K. In accord with our findings for isotope mixtures, 
we conclude that mean free path reduction by phonon- 
phonon scattering starts dominating the adverse effect of 
defects at high temperatures. 



B. Thermal conductivity of graphene 

Due to the present interest in graphene, we determined 
the influence of defects on its thermal conductivity and 
present our results in Fig. [2] Results for A in struc- 
turally perfect, free-standing infinite graphene monolay- 
ers with the isotopic composition l2, C x l2 Ci- x are shown 
in Figs. |2|a) and (b) for a selected set of compositions in 
the temperature range up to 600 K. Graphene monolay- 
ers were represented by a periodic array of rectangular 
180-atom unit cells. As in the case of nanotubes, isotopi- 
cally pure graphene has the highest thermal conductivity. 
Similar to nanotubes, the maximum value of A is reached 
near T^lOO K. 

We find the thermal conductivity to be slightly higher 
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FIG. 2. (Color online) Thermal conductivity A of perfect and defective graphene and graphene nanoribbons as a function of 
temperature T. (a) A in 13 C x 12 Ci- x graphene with varying isotopic composition in comparison to pure 12 C graphene with 
1% missing atoms forming divacancy defects, (b) Details of (a) on a reduced temperature scale, (c) A of 12 C graphene with 
divacancies, presented in (a) and (b), on an expanded A scale, (d) Depiction of the graphene unit cells containing different 
types of defects. The lines in (a)-(c) are guides to the eye. 



in defect-free graphene consisting of 12 C than of 13 C. In 
contrast to nanotubes, the reduction of the thermal con- 
ductivity in isotopic mixtures is much more pronounced 
in graphene. Whereas a 10% content of 13 C isotopic 
impurities reduced the thermal conductivity at 100 K 
by 30% in nanotubes, the corresponding 50% reduction 
in graphene is much larger. This is consistent with the 
fact that in defect-free systems, phonon-phonon scatter- 
ing limits the phonon mean free path more in nanotubes 
with a finite circumference than in graphene. 

Similar to nanotubes, the presence of a single diva- 
cancy per unit cell quenches the thermal conductivity by 
more than an order of magnitude, as seen in Figs. [2|a)- 
[5|c). Also in graphene, the reduction of the thermal con- 
ductivity due to isotopic and structural defects is mainly 
caused by the decrease in the phonon mean free path. 



For isotopically pure 12 C based systems, comparison 
between A of graphene in Fig.[2|a) and graphene nanorib- 
bons in Fig. [3] reveals that edges in nanoribbons quench 
thermal conductivity in a much more drastic way than 
sparse diatomic vacancies in an infinite graphene mono- 
layer. We observe reduction of A by a factor of 700 at 
T = 100 K, a factor of 30 at 300 K, and by an order of 
magnitude at 600K. 

Possibly unexpected at the first glance is our finding 
that placing one divacancy per 600-atom unit cell, cor- 
responding to 0.33% atomic vacancies, does not cause as 
drastic a reduction of the thermal conductivity as in the 
case of nanotubes and graphene. Whereas this effect is 
still significant at low temperatures, amounting to a 60% 
reduction at T = 100 K, it becomes negligibly small at 



C. Thermal conductivity of graphene nanorribons 

Finally, we compared thermal conductivity of graphene 
to that of 11.1 A wide armchair graphene nanoribbons, 
using the same 13 C x 12 Ci- x isotopic compositions as for 
graphene and nanotubes. The nanoribbons were repre- 
sented using periodic boundary conditions using 131.5 A 
long unit cells containing 600 C atoms. Our results for 
the thermal conductivity of nanoribbons are presented in 
Fig-i 

Our results for finite-width nanoribbons should be very 
relevant also for graphene formed by chemical vapor de- 
position (CVD). Grain boundaries in CVD graphene, 
viewed as lines of incorrectly coordinated carbon atoms, 
will scatter phonons and limit thermal conductivity in 
a very similar way as nanoribbon edges. Since the unit 
cell size representing realistic grain boundaries is pro- 
hibitively large for atomistic simulations, we only point 
out the analogy between nanoribbons and polycrystalline 
graphene. 



600 
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FIG. 3. (Color online) Temperature dependence of the 
thermal conductivity A of 13 C x 12 Ci- x structurally perfect 
graphene nanoribbons in comparison to pure 12 C graphene 
nanoribbons with 0.33% missing atoms arranged as divacan- 
cies. The lines are guides to the eye. 
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temperatures above 400 K. We conclude that the under- 
lying reduction of the phonon mean free path in narrow 
graphene nanoribbons is dominated by the presence of 
edges and that additional structural defects play only a 
minor role, especially at higher temperatures. 

Due to the dominating role of edges as scattering 
centers in nanoribbons, the role of isotopic impurities 
is much smaller in these systems than in nanotubes 
and graphene. Still, in agreement with our results for 
nanotubes and graphene, we find that isotopically pure 
nanoribbons based on 12 C conduct heat slightly better 
than those based on 13 C. 



D. High-temperature behavior 

Even though infinitely extended graphene appears to 
conduct heat better than nanotubes at low temperatures, 
the difference between the two systems becomes smaller 
at 400 K and above. As already mentioned above, the 
adverse effect of defects on the thermal conductivity of 
carbon nanostructures becomes less significant at higher 
temperatures, when uncorrelated atomic motion reduces 
the phonon mean free path even in defect-free systems. 
In practice, we find very similar thermal conductivities 
in isotopically pure systems and in isotope mixtures at 
very high temperatures T>600 K. At still higher tem- 
peratures approaching the melting point, when vacancy 
production occurs naturally, presence of additional struc- 
tural defects should play a negligible role as well. At 
those high temperatures, thermal conductivity of nan- 
otubes and graphene may drop close to that of nanorib- 
bons. 



IV. SUMMARY AND CONCLUSIONS 

In conclusion, we studied the adverse role of defects on 
the thermal conductivity of carbon nanotubes, graphene 
and graphene nanoribbons using non-equilibrium molec- 



ular dynamics simulations. We found that all defects, 
including divacancies, extended edges and isotopic im- 
purities reduce thermal conductivity significantly in all 
systems by introducing phonon scattering centers and 
thus decreasing the phonon mean free path. We rec- 
onciled results of former studies, which differed by up to 
an order of magnitude, by identifying limitations of var- 
ious computational approaches. We found that infinite, 
defect-free graphene should conduct heat better than any 
other carbon nanostructure at low temperatures. For 
temperatures T<400 K, isotopic impurities were found 
to quench the thermal conductivity of graphene more 
than that of carbon nanotubes. We found that even sub- 
percent concentrations of divacancies reduced the ther- 
mal conductivity of all nanocarbons more than much 
higher concentrations of isotopic impurities. For tem- 
peratures T<400 K, the adverse effect of divacancies was 
found to be more pronounced in graphene than in car- 
bon nanotubes. Finite-width graphene nanoribbons can 
be viewed as graphene with extended vacancies and thus 
have not only a much lower thermal conductivity, but 
also a lower susceptibility to the presence of additional 
defects than graphene. At high temperatures, when an- 
harmonicities in the force field reduce the phonon mean 
free path more than defects, we find that thermal con- 
ductivity decreases significantly and that differences be- 
tween particular nanocarbons become washed out to a 
large degree. 
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